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Abstract 

We develop a Bayesian "sum-of-trees" model where each tree is constrained by a 
regularization prior to be a weak learner, and fitting and inference are accomplished 
via an iterative Bayesian backfitting MCMC algorithm that generates samples from a 
posterior. Effectively, BART is a nonparametric Bayesian regression approach which 
uses dimensionally adaptive random basis elements. Motivated by ensemble methods in 
general, and boosting algorithms in particular, BART is defined by a statistical model: 
a prior and a likelihood. This approach enables full posterior inference including point 
and interval estimates of the unknown regression function as well as the marginal effects 
of potential predictors. By keeping track of predictor inclusion frequencies, BART can 
also be used for model free variable selection. BART's many features are illustrated 
with a bake-off against competing methods on 42 different data sets, with a simulation 
experiment and on a drug discovery classification problem. 
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1 Introduction 



We consider the fundamental problem of making inference about an unknown 
function / that predicts an output Y using a p dimensional vector of inputs 
x = (xx, . . . , Xp) when 

Y = f{x) + e, e~iV(0,a 2 ). (1) 

To do this, we consider modelling or at least approximating f(x) = E(Y \ x), the 
mean of Y given x, by a sum of m regression trees f(x) ~ h(x) = Y?jLi 9j{x) 
where each gi denotes a regression tree. Thus we approximate (1) by a sum-of- 
trees model 

Y = h(x) + e, e~N(0,a 2 ). (2) 

A sum-of-trees model is fundamentally an additive model with multivariate 
components. Compared to generalized additive models based on sums of low 
dimensional smoothers, these multivariate components can more naturally incor- 
porate interaction effects. And compared to a single tree model, the sum-of-trees 
can more easily incorporate additive effects. 

Various methods which combine a set of tree models, so called ensemble meth- 
ods, have attracted much attention. These include boosting (Freund & Schapire 
(1997), Friedman (2001)), random forests (Breiman 2001) and bagging (Breiman 
1996), each of which use different techniques to fit a linear combination of trees. 
Boosting fits a sequence of single trees, using each tree to fit data variation not ex- 
plained by earlier trees in the sequence. Bagging and random forests use data ran- 
domization and stochastic search to create a large number of independent trees, 
and then reduce prediction variance by averaging predictions across the trees. Yet 
another approach that results in a linear combination of trees is Bayesian model 
averaging applied to the posterior arising from a Bayesian single-tree model as 
in Chipman, George & McCulloch (1998) (hereafter CGM98), Denison, Mallick 
& Smith (1998) and Wu, Tjelmeland & West (2007). Such model averaging uses 
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posterior probabilities as weights for averaging the predictions from individual 
trees. 

In this paper we propose a Bayesian approach called BART (Bayesian Ad- 
ditive Regression Trees) which uses a sum of trees to model or approximate 
f(x) = E(Y | x). The essential idea is to elaborate the sum-of-trees model (2) 
by imposing a prior that regularizes the fit by keeping the individual tree effects 
small. In effect, the g/s become a dimensionally adaptive random basis of "weak 
learners", to borrow a phrase from the boosting literature. By weakening the gj 
effects, BART ends up with a sum of trees, each of which explains a small and 
different portion of /. Note that BART is not equivalent to posterior averaging 
of single tree fits of the entire function /. 

To fit the sum-of-trees model, BART uses a tailored version of Bayesian back- 
fitting MCMC (Hastie & Tibshirani 2000) that iteratively constructs and fits 
successive residuals. Although similar in spirit to the gradient boosting approach 
of Friedman (2001), BART differs in both how it weakens the individual trees by 
instead using a prior, and how it performs the iterative fitting by instead using 
Bayesian backfitting on a fixed number of trees. Conceptually, BART can be 
viewed as a Bayesian nonparametric approach that fits a parameter rich model 
using a strongly influential prior distribution. 

Inferences obtained from BART are based on successive iterations of the back- 
fitting algorithm which are effectively an MCMC sample from the induced pos- 
terior over the sum-of-trees model space. A single posterior mean estimate of 
f(x) = E(Y | x) at any input value x is obtained by a simple average of these 
successive sum-of-trees model draws evaluated at x. Further, pointwise uncer- 
tainty intervals for f(x) are easily obtained from the corresponding quantiles 
of the sample of draws. Point and interval estimates are similarly obtained for 
functionals of /, such as partial dependence functions which reveal the marginal 
effects of the x components. Finally, by keeping track of the relative frequency 
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with which x components appear in the sum-of-trees model iterations, BART 
can be used to identify which components are more important for explaining the 
variation of Y. Such variable selection information is model- free in the sense that 
it is not based on the usual assumption of an encompassing parametric model. 

Finally, to facilitate the use of the BART methods described in this paper, we 
have provided open-source software implementing BART as a stand-alone pack- 
age or with an interface to R, along with full documentation and examples. It is 
available at http : / /gsbwww . uchicago . edu/f ac/ robert . mcculloch/ research, 
and as the BayesTree library in R at http://cran.r-project.org/. 

The remainder of the paper is organized as follows. In Section 2, the BART 
model is outlined. This consists of the sum-of-trees model combined with a 
regularization prior. In Section 3, a Bayesian backfitting MCMC algorithm and 
methods for inference are described. In Section 4, we describe a probit extension 
of BART for classification of binary Y. In Section 5, examples, both simulated 
and real, are used to demonstrate the potential of BART. Section 6 describes 
extensions and a variety of recent developments and applications of BART based 
on an early version of this paper. Section 7 concludes with a discussion. 

2 The BART Model 

As described in the introduction, the BART model consists of two parts: a sum- 
of-trees model and a regularization prior on the parameters of that model. We 
describe each of these in detail in the following subsections. 

2.1 A Sum-of-Trees Model 

To elaborate the form of the sum-of-trees model (2), we begin by establishing 
notation for a single tree model. Let T denote a binary tree consisting of a 
set of interior node decision rules and a set of terminal nodes, and let M = 



1 



{[ii, H2, ■ ■ ■ , Hb} denote a set of parameter values associated with each of the b 
terminal nodes of T. The decision rules are binary splits of the predictor space 
of the form {x G A} vs {x ^ A} where A is a subset of the range of x. These are 
typically based on the single components of x = (x±, . . . , x p ) and are of the form 
{xi < c} vs {xi > c} for continuous Xj. Each x value is associated with a single 
terminal node of T by the sequence of decision rules from top to bottom, and is 
then assigned the fii value associated with this terminal node. For a given T and 
M, we use g(x; T, M) to denote the function which assigns a ^ G M to x. Thus, 

Y = g(x;T,M) + e, e ~ N(0, a 2 ) (3) 

is a single tree model of the form considered by CGM98. Under (3), the condi- 
tional mean of Y given x, E(Y \x) equals the terminal node parameter /ij assigned 
by g(x;T, M). 

With this notation, the sum-of-trees model (2) can be more explicitly ex- 
pressed as 

Y= ff>(x;T,-,M,)j+e, e~iV(0,a 2 ), (4) 

where for each binary regression tree Tj and its associated terminal node pa- 
rameters Mj, g(x;Tj,Mj) is the function which assigns fcj G Mj to x. Under 
(4), E(Y | x) equals the sum of all the terminal node /%'s assigned to x by the 
g(x;Tj, Mj)'s. When the number of trees m > 1, each here is merely a part 
of E(Y | x), unlike the single tree model (3). Furthermore, each such /i^ will 
represent a main effect when g(x;Tj, Mj) depends on only one component of x 
(i.e., a single variable), and will represent an interaction effect when g(x; Tj, Mj) 
depends on more than one component of x (i.e., more than one variable). Thus, 
the sum-of-trees model can incorporate both main effects and interaction effects. 
And because (4) may be based on trees of varying sizes, the interaction effects 
may be of varying orders. In the special case where every terminal node assign- 
ment depends on just a single component of x, the sum-of-trees model reduces to 
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a simple additive function, a sum of step functions of the individual components 
of x. 

With a large number of trees, a sum-of-trees model gains increased representa- 
tion flexibility which, as we'll see, endows BART with excellent predictive capabil- 
ities. This representational flexibility is obtained by rapidly increasing the num- 
ber of parameters. Indeed, for fixed m, each sum-of-trees model (4) is determined 
by (Ti, Mi), . . . , (T m , M m ) and a, which includes all the bottom node parameters 
as well as the tree structures and decision rules. Further, the representational flex- 
ibility of each individual tree leads to substantial redundancy across the tree com- 
ponents. Indeed, one can regard {g(x; Ti, Mi), . . . , g{x; T m , M m )} as an "over- 
complete basis" in the sense that many different choices of (Ti, Mi), . . . , (T m , M m ) 
can lead to an identical function J2]Li g{ x ', Tj, Mj). 

2.2 A Regular izat ion Prior 

We complete the BART model specification by imposing a prior over all the pa- 
rameters of the sum-of-trees model, namely (Ti, Mi), . . . , (T m , M m ) and a. As 
discussed below, we advocate specifications of this prior that effectively regularize 
the fit by keeping the individual tree effects from being unduly influential. With- 
out such a regularizing influence, large tree components would overwhelm the rich 
structure of (4), thereby limiting the advantages of the additive representation 
both in terms of function approximation and computation. 

To facilitate the easy implementation of BART in practice, we recommend 
automatic default specifications below which appear to be remarkably effective, 
as demonstrated in the many examples of Section 5. Basically we proceed by 
first reducing the prior formulation problem to the specification of just a few 
interpretable hyperparameters which govern priors on Tj, Mj and a. Our rec- 
ommended defaults are then obtained by using the observed variation in y to 
gauge reasonable hyperparameter values when external subjective information 
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is unavailable. Alternatively one can use the computationally more demanding 
approach of hyperparameter selection via cross validation. 



2.2.1 Prior Independence and Symmetry 



Specification of our regularization prior is vastly simplified by restricting atten- 
tion to priors for which 



p{(T 1 ,M 1 ),...,{T m ,M m ),a) 



p(a) 



]Jp( M j I T j)p( T j 



p(a) 



and 



p(M j \T J ) = Y[p( f i tl \T 3 ), 



(5) 



(6) 



where Hij G Mj. Under such priors, the tree components (Tj, Mj) are independent 
of each other and of a, and the terminal node parameters of every tree are 
independent. 

The independence restrictions above simplify the prior choice problem to the 
specification of prior forms for just p(Tj),p(fiij \ Tj) and p(cr), and to simplify 
it further we consider identical forms for all p{Tj) and for all | Tj). As 
described in the ensuing subsections, we use the same prior forms for these as 
those proposed by CGM98 for Bayesian CART. In addition to their valuable 
computational benefits, these forms are controlled by just a few interpretable 
hyperparameters which can be calibrated using the data to yield effective default 
specifications for regularization of the sum-of-trees model. However, as will be 
seen, considerations for the choice of these hyperparameter values for BART are 
markedly different than those for Bayesian CART. 
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2.2.2 The 7} Prior 

For p(Tj), the form recommended by CGM98 is easy to specify and dovetails 
nicely with calculations for the backfitting MCMC algorithm described later in 
Section 3.1. It is specified by three aspects: (i) the probability that a node at 
depth d is nonterminal, given by 

a{l + d)-^ a E (0, 1), (3 G [0, oo), (7) 

(ii) the distribution on the splitting variable assignments at each interior node, 
and (iii) the distribution on the splitting rule assignment in each interior node, 
conditional on splitting variable. For (ii) and (iii) we use the simple defaults 
used by CGM98, namely the uniform prior on available variables for (ii) and the 
uniform prior on the discrete set of available splitting values for (iii). 

In a single tree model, (i.e. m — 1), a tree with many terminal nodes may 
be needed to model complicated structure. However, for a sum-of-trees model, 
especially with m large, we want the regularization prior to keep the individual 
tree components small. In our examples in Section 5, we do so by using a = .95 
and f3 — 2 in (7). With this choice, trees with 1, 2, 3, 4, and > 5 terminal nodes 
receive prior probability of 0.05, 0.55, 0.28, 0.09, and 0.03, respectively. Note that 
even with this prior, which puts most probability on tree sizes of 2 or 3, trees 
with many terminal nodes can be grown if the data demands it. For example, in 
one of our simulated examples with this prior, we observed considerable posterior 
probability on trees of size 17 when we set m = 1. 

2.2.3 The ^ | Tj Prior 

For p(fiij | Tj), we use the conjugate normal distribution iV^/i^cx^) which offers 
tremendous computational benefits because fiij can be margined out. To guide 
the specification of the hyperparameters /i M and <r M , note that E(Y \ x) is the sum 
of m /iy 's under the sum-of-trees model, and because the /x^'s are apriori iid, the 
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induced prior on E(Y | x) is N(m fi^,ma^). Note also that it is highly probable 
that E(Y | x) is between y m i n and y ma x, the observed minimum and maximum of 
Y in the data. The essence of our strategy is then to choose /i M and <r M so that 
N(m fi^jma^) assigns substantial probability to the interval (y m in,ymax)- This 
can be conveniently done by choosing /i M and <r M so that mfi^ — k y/ma^ = y m i n 
and mix^ + k ^Jma l^ = y m ax for some preselected value of k. For example, k = 2 
would yield a 95% prior probability that E(Y \ x) is in the interval (y m i n ,ymax)- 
The strategy above uses an aspect of the observed data, namely y m in and y ma x, 
to try to ensure that the implicit prior for E(Y | x) is in the right "ballpark". 
That is to say, we want it to assign substantial probability to the entire region of 
plausible values of E(Y | x) while avoiding overconcentration and overdispersion. 
We have found that, as long as this goal is met, BART is very robust to our 
exact specification. Such a data-informed prior approach is especially useful in 
our problem, where reliable subjective information about E(Y | x) is likely to be 
unavailable. 

For convenience, we implement our specification strategy by first shifting and 
rescaling Y so that the observed transformed y values range from y m in = —0.5 
to y m ax — 0.5, and then treating this transformed Y as our dependent variable. 
We then simply center the prior for fiy at zero /i M = and choose so that 
ky/mcr^ = 0.5 for a suitable value of k, yielding 



This prior has the effect of shrinking the tree parameters /i^ towards zero, 
limiting the effect of the individual tree components of (4) by keeping them small. 
Note that as k and/or the number of trees m is increased, this prior will become 
tighter and apply greater shrinkage to the fii/s. Prior shrinkage on /ij/s is the 
counterpart of the shrinkage parameter in Friedman's (2001) gradient boosting 
algorithm. The prior standard deviation of /i^ here and the gradient boosting 
shrinkage parameter there, both serve to "weaken" the individual trees so that 




8 



9 



each is constrained to play a smaller role in the overall fit. For the choice of 
k, we have found that values of k between 1 and 3 yield good results, and we 
recommend k = 2 as an automatic default choice. Alternatively the value of k 
may be chosen by cross-validation. 

Although the calibration of this prior is based on a simple linear transforma- 
tion of Y, it should be noted that there is no need to transform the predictor 
variables. This is a consequence of the fact that the tree splitting rules are in- 
variant to monotone transformations of the x components. The simplicity of our 
prior for fiy is an appealing feature of BART. In contrast, methods like neural 
nets that use linear combinations of predictors require standardization choices 
for each predictor. 

2.2.4 The a Prior 

For p(cr), we also use a conjugate prior, here the inverse chi-square distribution 
a 2 ~ v X/xt- To guide the specification of the hyperparameters v and A, we again 
use a data-informed prior approach, in this case to assign substantial probability 
to the entire region of plausible values of a while avoiding overconcentration and 
overdispersion. Essentially, we calibrate v and A for this purpose using a "rough 
data-based overestimate" a of a. Two natural choices for a are 1) the "naive" 
specification, in which we take a to be the sample standard deviation of Y, or 
2) the "linear model" specification, in which we take a as the residual standard 
deviation from a least squares linear regression of Y on the original X's. We then 
pick a value of v between 3 and 10 to get an appropriate shape, and a value of A 
so that the gth quantile of the prior on a is located at a, that is P(a < a) = q. 
We consider values of q such as 0.75, 0.90 or 0.99 to center the distribution below 

(7. 

Figure 1 illustrates priors corresponding to three (z/, q) settings when the rough 
overestimate is o = 2. We refer to these three settings, (V, q) = (10,0.75), 
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Figure 1: Three priors on a when a = 2. 

(3,0.90), (3,0.99), as conservative, default and aggressive, respectively. The 
prior mode moves towards smaller a values as q is increased. We recommend 
against choosing v < 3 because it seems to concentrate too much mass on very 
small a values, which leads to overfitting. In our examples, we have found these 
three settings work very well and yield similar results. For automatic use, we 
recommend the default setting (u, q) = (3, 0.90) which tends to avoid extremes. 
Alternatively, the values of (u, q) may be chosen by cross-validation. 

2.2.5 The Choice of m 

A major difference between BART and boosting methods is that for a fixed 
number of trees m, BART uses an iterative backfitting algorithm (described in 
Section 3.1) to cycle over and over through the m trees. If BART is to be used 
for estimating f(x) or predicting Y, it might be reasonable to treat m as an 
unknown parameter, putting a prior on m and proceeding with a fully Bayes im- 
plementation of BART. Another reasonable strategy might be to select a "best" 
value for m by cross-validation. However, both of these strategies substantially 
increase computational requirements. 



- - - conservative: df=10, quantile=.75 

default: df=3, quantile=.9 

- aggressive: df=3, quantile=.99 
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To avoid the computational costs of these strategies, we have found it fast 
and expedient for estimation and prediction to begin with a default of m = 200, 
and then perhaps to check if one or two other choices makes any difference. Our 
experience has been that as m is increased, starting with m = 1, the predictive 
performance of BART improves dramatically until at some point it levels off and 
then begins to very slowly degrade for large values of m. Thus, for prediction, it 
seems only important to avoid choosing m too small. As will be seen in Section 
5, BART yielded excellent predictive performance on a wide variety of examples 
with the simple default m = 200. 

However, as mentioned in Section 1, BART can also be used for variable se- 
lection by simply selecting those variables that appear most often in the fitted 
sum-of-trees models. Interestingly, this strategy does not seem to work so well 
when m is large, because the sum-of-trees model incorporates many irrelevant 
predictions in addition to the relevant ones, a consequence of the rich flexibility 
granted by the redundancy of having so many trees. However, as m is decreased 
and that redundancy is diminished, BART tends to heavily favor relevant pre- 
dictors for its fit. In a sense, when m is small the predictors compete with each 
other to improve fit so that those that appear more often are likely to be more 
useful for prediction. As discussed in Section 3.2 and illustrated in Section 5, a 
useful implementation of this strategy is to use BART as an exploratory tool by 
observing what happens to component inclusion frequencies as m is decreased. 

3 Extracting Information from the Posterior 

3.1 A Bayesian Backfitting MCMC Algorithm 

Given the observed data y, our Bayesian setup induces a posterior distribution 

p((Ti, Mi), . . . , (T m , M m ), a\ y) (9) 
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on all the unknowns that determine a sum-of-trees model (4). Although the 
sheer size of the parameter space precludes exhaustive calculation, the following 
backfitting MCMC algorithm can be used to sample from this posterior. 

At a general level, our algorithm is a Gibbs sampler. For notational conve- 
nience, let Tq-) be the set of all trees in the sum except Tj, and similarly define 
M(j). Thus Tq-) will be a set of m — 1 trees, and Mq) the associated terminal 
node parameters. The Gibbs sampler here entails m successive draws of (Tj, Mj) 
conditionally on (T(j), Mq), a): 

(T^Mj^T^M^^y, (10) 

j = 1, . . . , m, followed by a draw of a from the full conditional: 

a\T l ,...T m ,M 1 ,...,M m ,y. (11) 

Hastie & Tibshirani (2000) considered a similar application of the Gibbs sampler 
for posterior sampling for additive and generalized additive models with a fixed, 
and showed how it was a stochastic generalization of the backfitting algorithm 
for such models. For this reason, we refer to our algorithm as backfitting MCMC. 

The draw of a in (11) is simply a draw from an inverse gamma distribution 
and so can be easily obtained by routine methods. More challenging is how to 
implement the m draws of (Tj,Mj) in (10). This can be done by taking advan- 
tage of the following reductions. First, observe that the conditional distribution 
p(Tj, Mj\T {j) , M (j) , a, y) depends on (T {j) , M (j) , y) only through 

Rj = y-^2g(x;T k ,M k ), (12) 

the n— vector of partial residuals based on a fit that excludes the j'th tree. Thus, 
the m draws of (Tj, Mj) given (Tq), M^, a, y) in (10) are equivalent to m draws 
from 

(Tj,Mj)\Rj,a, (13) 
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j = 1, . . . ,m. 

Now (13) is formally equivalent to the posterior of the single tree model Rj = 
g(x; Tj, Mj) + e where Rj plays the role of the data y. Because we have used a 
conjugate prior for Mj, 

p(T j \R j ,a)<xp(T j ) J p{Rj\Mj,Tj,a)p{Mj\Tj,a)dMj (14) 

can be obtained in closed form up to a norming constant. This allows us to carry 
out each draw from (13) in two successive steps as 

Tj\Rj,a (15) 

Mj\Tj,Rj,a. (16) 

The draw of Tj in (15), although somewhat elaborate, can be obtained using 
the Metropolis-Hastings (MH) algorithm of CGM98. This algorithm proposes a 
new tree based on the current tree using one of four moves. The moves and their 
associated proposal probabilities are: growing a terminal node (0.25), pruning a 
pair of terminal nodes (0.25), changing a non-terminal rule (0.40), and swapping 
a rule between parent and child (0.10). Although the grow and prune moves 
change the number of terminal nodes, by integrating out Mj in (14), we avoid 
the complexities associated with reversible jumps between continuous spaces of 
varying dimensions (Green 1995). 

Finally, the draw of Mj in (16) is simply a set of independent draws of the 
terminal node fcj's from a normal distribution. The draw of Mj enables the cal- 
culation of the subsequent residual Rj+i which is critical for the next draw of Tj. 
Fortunately, there is again no need for a complex reversible jump implementation. 

We initialize the chain with m simple single node trees, and then iterations 
are repeated until satisfactory convergence is obtained. At each iteration, each 
tree may increase or decrease the number of terminal nodes by one, or change 
one or two decision rules. Each \i will change (or cease to exist or be born), and 
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cr will change. It is not uncommon for a tree to grow large and then subsequently 
collapse back down to a single node as the algorithm iterates. The sum-of- 
trees model, with its abundance of unidentified parameters, allows for "fit" to 
be freely reallocated from one tree to another. Because each move makes only 
small incremental changes to the fit, we can imagine the algorithm as analogous 
to sculpting a complex figure by adding and subtracting small dabs of clay. 

Compared to the single tree model MCMC approach of CGM98, our backfit- 
ting MCMC algorithm mixes dramatically better. When only single tree models 
are considered, the MCMC algorithm tends to quickly gravitate towards a single 
large tree and then gets stuck in a local neighborhood of that tree. In sharp 
contrast, we have found that restarts of the backfitting MCMC algorithm give 
remarkably similar results even in difficult problems. Consequently, we run one 
long chain with BART rather than multiple starts. Although mixing does not 
appear to be an issue, the recently proposed modifications of Wu et al. (2007) 
might well provide additional benefits. 

3.2 Posterior Inference Statistics 

The backfitting algorithm described in the previous section is ergodic, generating 
a sequence of draws of (T 1; Ml), . . . , (T m , M m ), a which is converging (in distri- 
bution) to the posterior p((Ti, Mi), . . . , (T m , M m ), a\ y). The induced sequence 
of sum-of-trees functions 

m 

/*(•) = £<?(- ; t;,m;), (i?) 

3=1 

for the sequence of draws (T x *, M^), . . . , (T^, M^J, is thus converging to p(f \ y), 
the posterior distribution on the "true" /(•). Thus, by running the algorithm long 
enough after a suitable burn-in period, the sequence of /* draws, say /*,..., /£, 
may be regarded as an approximate, dependent sample from p(f | y). Bayesian 
inferential quantities of interest can then be approximated with this sample as 
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indicated below. Although the number of iterations needed for reliable inferences 
will of course depend on the particular application, our experience with the ex- 
amples in Section 5 suggests that the number of iterations required is relatively 
modest. 

To estimate f(x) or predict Y at a particular x, in-sample or out-of-sample, 
a natural choice is the average of the after burn-in sample fl , . . . , /£, 

^E/M (is) 
1X k=i 

which approximates the posterior mean E(f(x) \ y). Another good choice would 
be the median of f*(x), . . . , f^(x) which approximates the posterior median 
of f{x). Posterior uncertainty about f(x) may be gauged by the variation of 
fi(x), . . . , fx(x). For example, a natural and convenient (1 — a)% posterior 
interval for f(x) is obtained as the interval between the upper and lower a/2 
quantiles of f*(x), . . . , fx{x). As will be seen, these uncertainty intervals behave 
sensibly, for example by widening at x values far from the data. 

It is also straightforward to use fi(x), . . . , fx{x) to estimate other functionals 
of /. For example, a functional of particular interest is the partial dependence 
function, (Friedman 2001), which summarizes the marginal effect of one (or more) 
predictors on the response. More precisely, letting f(x) = f(x s ,x c ) where x has 
been partitioned into the predictors of interest, x s and the complement x c = x\x s , 
the partial dependence function is defined as 

1 n 

f( x s) = - f( x *' x ^), (19) 

where Xj C is the ith observation of x c in the data. Note that (x s ,Xi C ) will not 
generally be one of the observed data points. A draw from the induced BART 
posterior p(f(x s ) \ y) at any value of x s is obtained by simply computing f%(x s ) = 
fk( x si x ic)- The average of f*(x s ), ■ ■ ■ , h* K (x s ) then yields an estimate of 
f(x s ), and the upper and lower a/2 quantiles provide endpoints of (1 — a)% 
posterior intervals for f(x s ). 



16 



As discussed in Section 2.2.5, BART can also be used as an exploratory tool 
for variable selection. This can be accomplished by observing what happens to 
the x component frequencies in a sequence of MCMC samples f*, ■ ■ ■ , fx as the 
number of trees m is set smaller and smaller. More precisely, for each simulated 
sum-of-trees model f£, let be the proportion of all splitting rules that use the 
ith component of x. Then 

1 K 

Vi = -J2zik (20) 
n k=i 

is the average use per splitting rule for the ith component of x. As m is set 
smaller and smaller, the sum-of-trees models tend to more strongly favor inclusion 
of those x components which improve prediction of y and exclusion of those x 
components that are unrelated to y. In effect, smaller m seems to create a 
bottleneck that forces the x components to compete for entry into the sum-of- 
trees model. As is illustrated in Section 5, the x components with the larger fj's 
will then be those that provide the most information for predicting y. 



4 Extending BART for Classification 

Our development of BART up to this point pertains to setups where the output 
of interest Y is a continuous variable. However, for binary Y (— or 1), it is 
straightforward to extend BART to the probit model setup for classification 

p(x) = P[Y = 1 1 x] = &[G(x)] (21) 

where 

m 

G(x) = Y:9(x;T j ,M j ) (22) 
i=i 

and $[■] is the standard normal cdf. For this extension of BART, we need to 
impose a regularization prior on G(x) and to implement a Bayesian backfitting 
algorithm for posterior computation. Fortunately, these are obtained with only 
minor modifications of the methods in Sections 2 and 3. 
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As opposed to (4), the model (21) does not involve a and so we need only put 
a prior on only (7\, Mi), . . . , (T m , M m ). Proceeding exactly as in Section 2.2.1, 
we consider a prior of the form 



p((T 1; Mi), . . . , (T m , M m )) = [] p(Tj) UpifHj | Tj) 



(23) 



3 L * 



where each tree prior p(Tj) is the choice recommended in Section 2.2.2. For the 
choice of | Tj) here, we consider the case where the interval ($[—3.0], $[3.0]) 
contains most of the p(x) values of interest, a case which will often be of practical 
relevance. Proceeding similarly to the motivation of (8) in Section 2.2.3, we would 
then recommend the choice 



where k is such that G{x) will with high probability be in the interval (—3.0, 3.0). 
Just as for (8), this prior has the effect of shrinking the tree parameters pij 
towards zero, limiting the effect of the individual tree components of G(x). As 
k and/or the number of trees m is increased, this prior will become tighter and 
apply greater shrinkage to the /ij/s. For the choice of k, we have found that 
values of k between 1 and 3 yield good results, and we recommend k = 2 as 
an automatic default choice. Alternatively the value of k may be chosen by 
cross-validation. 

By shrinking G(x) towards 0, the prior (24) has the effect of shrinking p(x) = 
Q[G(x)} towards 0.5. If it is of interest to shrink towards a value po other than 
0.5, one can simply replace G(x) by G c = G(x) + c in (21) where c = $~ 1 [j>o]. 
Note also that if an interval other than ($[—3.0], $[3.0]) is of interest for p(x), 
suitable modification of (24) is straightforward. 

Turning to posterior calculation, the essential features of the backfitting al- 
gorithm in Section 3.1 can be implemented by using the augmentation idea of 
Albert & Chib (1993). The key idea is to recast the model (21) by introducing 



fj4j ~ A r (0, £7^) where = 3.0/ky/m 



(24) 
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latent variables Z 1; . . . , Z n iid ~ N(G(x), 1) such that = 1 if Z± > and Y"j = 
if Zi < 0. Note that under this formulation, | = 1] ~ max{iV(g(£), 1),0} 
and Zi \ [yi = 0] ~ min{N(g(x), 1), 0}. Incorporating simulation of the latent 
Zi values into the backfitting algorithm, the Gibbs sampler iterations here en- 
tail n successive draws of Zi\y^ i = 1, . . . , n followed by m successive draws 
of (Tj, Mj)\T(j), My), Zx, z n , j = l,...,m, as spelled out in Section 3.1. The 
induced sequence of sum-of-trees functions 



p*(.)=* 



3=1 



(25) 



for the sequence of draws (T*, M±), . . . , (T^, M^), is thus converging to the poste- 
rior distribution on the "true" p(-). After a suitable burn-in period, the sequence 
of g* draws, say . . . , g* K , may be regarded as an approximate, dependent sam- 
ple from this posterior which can be used to draw inference about p(-) in the 
same way that /*,...,/# was used in Section 3.2 to draw inference about /(•). 



5 Applications 

In this section we demonstrate the application of BART on several examples. We 
begin in Section 5.1 with a predictive cross-validation performance comparison 
of BART with competing methods on 42 different real data sets. We next, in 
Section 5.2, evaluate and illustrate BART's capabilities on simulated data used 
by Friedman (1991). Finally, in Section 5.3 we apply the BART probit model to a 
drug discovery classification problem. All of the BART calculations throughout 
this section can be reproduced with the BayesTree library at http : //cran . 
r-project . org/. 
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5.1 Predictive Comparisons on 42 Data Sets 

Our first illustration is a "bake-off", a predictive performance comparison of 
BART with competing methods on 42 different real data sets. These data sets 
are a subset of 52 sets considered by Kim, Loh, Shih & Chaudhuri (2007). Ten 
data sets were excluded either because Random Forests was unable to use over 
32 categorical predictors, or because a single train/test split was used in the 
original paper. All data sets correspond to regression setups with between 3 and 
28 numeric predictors and to 6 categorical predictors. Categorical predictors 
were converted into 0/1 indicator variables corresponding to each level. Sample 
sizes vary from 96 to 6806 observations. In each of the 42 data sets, the response 
was minimally preprocessed, applying a log or square root transformation if this 
made the histogram of observed responses more bell-shaped. In about half the 
cases, a log transform was used to reduce a right tail. In one case (Fishery) a 
square root transform was most appropriate. 

For each of the 42 data sets, we created 20 independent train/test splits by 
randomly selecting 5/6 of the data as a training set and the the remaining 1/6 
as a test set. Thus, 42 x 20 = 840 test/train splits were created. Based on each 
training set, each method was then used to predict the corresponding test set 
and evaluated on the basis of its predictive RMSE. 

We considered two versions of BART: BART-cv where the prior hyperparam- 
eters (u, q, k, m) were treated as operational parameters to be tuned via cross- 
validation, and BART-default where we set (v, q, k, m) to the defaults (3, 0.90, 2, 200). 
For both BART-cv and BART-default, all specifications of the quantile q were 
made relative to the least squares linear regression estimate a, and the number 
of burn-in steps and MCMC iterations used were determined by inspection of a 
single long run. Typically 200 burn- in steps and 1000 iterations were used. For 
BART prediction at each x, we used the posterior mean estimates given by (18). 

As competitors we considered linear regression with LI regularization (the 
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Lasso) (Efron, Hastie, Johnstone & Tibshirani 2004) and three black-box mod- 
els: gradient boosting (Friedman (2001), implemented as gbm in R by Ridgeway 
(2004)), random forests (Breiman 2001), and neural networks with one layer 
of hidden units, (implemented as nnet in R by Venables & Ripley (2002)). 
These competitors were chosen because, like BART, they are black box predic- 
tors. Trees, Bayesian CART (CGM98), and Bayesian treed regression (Chipman, 
George & McCulloch 2002) models were not considered, since they tend to sac- 
rifice predictive performance for interpretability. 

With the exception of BART-default (which requires no tuning), the opera- 
tional parameters of every method were chosen via 5-fold cross-validation within 
each training set. The parameters considered and potential levels are given in 
Table 1. In particular, for BART-cv, we considered 

• three settings (3,0.90) (default), (3,0.99) (aggressive) and (10,0.75) (conser- 
vative) as shown in Figure 1 for the a prior hyperparameters (f, <?), 

• three values k = 2, 3, 5 reflecting moderate to heavy shrinkage for the \x 
prior hyperparameter, and 

• two values m = 50, 200 for the number of trees, 

a total of 3*3*2 = 18 potential choices for (z/, q, k, m). 

All the levels in Table 1 were chosen with a sufficiently wide range so that 
the selected value was not at an extreme of the candidate values in most prob- 
lems. Neural networks are the only model whose operational parameters need 
additional explanation. In that case, the number of hidden units was chosen 
in terms of the implied number of weights, rather than the number of units. 
This design choice was made because of the widely varying number of predictors 
across problems, which directly impacts the number of weights. A number of 
hidden units were chosen so that there was a total of roughly u weights, with 
u = 50, 100, 200, 500 or 800. In all cases, the number of hidden units was further 
constrained to fall between 3 and 30. For example, with 20 predictors we used 
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Method 


Parameter 


Values considered 


BART-cv 


Sigma prior: (i>, q) combinations 


(3,0.90), (3,0.99), (10,0.75) 




# trees m 


50, 200 




fi Prior: k value for cr M 


2, 3, 5 


Lasso 


shrinkage (in range 0-1) 


0.1, 0.2, 1.0 


Gradient 


# of trees 


50, 100, 200 


Boosting 


Shrinkage (multiplier of each tree 
added) 


0.01, 0.05, 0.10, 0.25 




Max depth permitted for each tree 


1,2, 3,4 


Neural 


# hidden units 


see text 


Nets 


Weight decay 


.0001, .001, .01, .1, 1, 2, 3 


Random 


# of trees 


500 


Forests 


% variables sampled to grow each 
node 


10, 25, 50, 100 



Table 1: Operational parameters for the various competing models. 



3, 8 and 21 as candidate values for the number of hidden units. 

To facilitate performance comparisons across data sets, we considered relative 
RMSE (RRMSE) which we defined as the RMSE divided by the minimum RMSE 
obtained by any method for each test/train split. Thus a method obtained an 
RRMSE of 1.0 when that method had the minimum RMSE on that split. As 
opposed to the RMSE, the RRMSE provides meaningful comparisons across data 
sets because of its invariance to location and scale transformations of the response 
variables. Boxplots of the 840 test /train split RRMSE values for each method 
are shown in Figure 2, and the median RRMSE values (the center of each box in 
Figure 2) are given in Table 2. (The Lasso was left off the boxplots because its 
many large RRMSE values visually overwhelmed the other comparisons). 

It is clear from the distribution of RRMSE values in Figure 2 that BART-cv 
tended to more often obtain smaller RMSE than any of its competitors. Per- 
haps even more interesting is that the overall performance of BART-default was 
arguably second best. This is especially impressive since neural nets, random 
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Figure 2: Boxplots of the RRMSE values for each method across the 840 test /train splits. 
Percentage RRMSE values larger than 1.5 for each method (and not plotted) were: Random 
forests 15.6%, Neural net 8.6%, Boosting 12.7%, BART-cv 9.0% and BART-default 11.2%. 
The Lasso (not plotted because of too many large RRMSE values), had 29.6% greater than 
1.5. 



Method Median 

Lasso 1.195 

Boosting 1.067 

BART-default 1.055 

Neural Net 1.055 

Random Forests 1.053 

BART-cv 1.039 



Table 2: Median RRMSE values for each method across the 840 test/train splits. 
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forests and gradient boosting all relied here on cross validation for control param- 
eter tuning. By avoiding the need for parameter specification, BART-default is 
vastly easier and faster to use. For example, a single implementation of BART-cv 
here requires selection among the 18 possible hyperparameter values with 5 fold 
cv, followed by fitting the best model, for a total of 18*5 + 1 = 91 applications 
of BART. Ready for easy "off the shelf" use, computationally inexpensive and 
coherent from a Bayesian point of view (cross validation is not), BART-default 
is in a sense the real winner in this experiment. It may also be of interest to 
note that, as opposed to BART-cv, posterior interval estimates based on BART- 
default retain their Bayesian validity. 

5.2 Friedman's Five Dimensional Test Function 

We next proceed to illustrate various features of BART on simulated data where 
we can gauge its performance against the true underlying signal. For this purpose, 
we constructed data by simulating values of x — (xi, x 2 , . . . , x p ) where 

xi, x 2 , ■ ■ ■ , x p iid ~ Uniform(0, 1), (26) 

and y given x where 

y = f( x ) + e = 10 sin(7rxix 2 ) + 20(x 3 - .5) 2 + 10x 4 + 5x 5 + e (27) 

where e ~ N(0, 1). Because y only depends onx\, . . . , x$, the predictors xq, . . . , x p 
are irrelevant. These added variables together with the interactions and nonlin- 
earities make it more challenging to find f(x) by standard parametric methods. 
Friedman (1991) used this setup with p = 10 to illustrate the potential of multi- 
variate adaptive regression splines (MARS). 

In Section 5.2.1, we illustrate various basic features of BART. We illustrate 
point and interval estimation of f(x), model-free variable selection and estima- 
tion of partial dependence functions. We see that the BART MCMC burns-in 
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quickly and mixes well. We illustrate BART's very robust performance with 
respect to hyperparameter changes. In Section 5.2.2, we show that BART out- 
performs competing methods including gradient boosting and neural nets. In 
Section 5.2.3, we elaborate this application to show BART's effectiveness at de- 
tecting low dimensional structure within high dimensional data. 

5.2.1 A Simple Application of BART 

We begin by illustrating the basic features of BART on a single simulated data 
set of the Friedman function (26) and (27) with p = 10 x's and n = 100 obser- 
vations. For simplicity, we applied BART with the default setting (u, q, k, m) = 
(3, 0.90, 2, 200) described in Section 2.2. Using the backfitting MCMC algorithm, 
we generated 5000 MCMC draws of /* as in (17) from the posterior after skipping 
1000 burn-in iterations. 

To begin with, for each value of x, we obtained posterior mean estimates 
f(x) of f(x) by averaging the 5000 f*(x) values as in (18). Endpoints of 90% 
posterior intervals for each f(x) were obtained as the 5% and 95% quantiles of the 
/* values. Figure 3(a) plots f(x) against f(x) for the n = 100 in-sample values 
of x from (26) which were used to generate the y values using (27). Vertical lines 
indicate the 90% posterior intervals for the /(x)'s. Figure 3(b) is the analogous 
plot at 100 randomly selected out-of-sample x values. We see that in-sample 
the f(x) values correlate very well with the true f(x) values and the intervals 
tend to cover the true values. Out-of sample, there is a slight degradation of the 
correlation and wider intervals indicating greater uncertainty about f(x) at new 
x values. 

Although we do not expect the 90% posterior intervals to exhibit 90% frequen- 
tist coverage, it may be of interest to note that 89% and 96% of the intervals in 
Figures 3(a) and (b) covered the true f(x) value, respectively. In fact, in over 200 
independent replicates of this example we found average coverage rates of 87% 
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(a) (b) (c) 




5 10 15 20 25 30 5 10 15 20 25 30 1000 3000 5000 

in-sample f(x) out-of-sample f(x) mcmc iteration 

Figure 3: Inference about Friedman's f(x) in p = 10 dimensions. 

(in-sample) and 93% (out-of-sample). In real data settings where / is unknown, 
bootstrap and/or cross-validation methods might be helpful to get similar cali- 
brations of frequentist coverage. In any case, however, it should be noted that 
for extreme x values, the prior may exert more shrinkage towards leading to 
lower coverage frequencies. 

The lower sequence in Figure 3(c) is the sequence of a draws over the entire 
1000 burn-in plus 5000 iterations (plotted with *). The horizontal line is drawn at 
the true value a — 1. The Markov chain here appears to reach equilibrium quickly, 
and although there is autocorrelation, the draws of a nicely wander around the 
true value a — 1 suggesting that we have fit but not overfit. To further highlight 
the deficiencies of a single tree model, the upper sequence (plotted with •) in 
Figure 3(c) is a sequence of a draws when m = 1, a single tree model, is used. 
The sequence seems to take longer to reach equilibrium and remains substantially 
above the true value a = 1. Evidently a single tree is inadequate to fit this data. 

Moving beyond estimation and inference about the values of f(x), BART 
estimates of the partial dependence functions f(xi) in (19) reveal the marginal 
effects of the individual x^s on y. Figure 4 shows the plots of point and inter- 
val estimates of the partial dependence functions for X\, . . . ,Xio from the 5000 
MCMC samples of /*. The nonzero marginal effects of x±, . . . ,x§ and the zero 
marginal effects of xq, . . . , x\o seem to be completely consistent with the form of 
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/ which of course would be unknown in practice. 

As described in Section 3.2, BART can also be used as a screening method 
for model-free variable selection information by keeping track of the average use 
per splitting rule v$ in (20) for each Xi over MCMC samples of /* for various 
choices of m. Figure 5 plots these i^'s for x±, . . . , X\o for m = 10, 20, 50, 100, 200. 
Notice that as the number of trees m is made smaller, the fitted sum-of-trees 
model tends to incorporate only those x variables, namely x\, . . . ,x$, that are 
needed to explain the variation of y. Thus, if one did not know the form of / 
here, selecting those variables which appear most often in the fitted sum-of-tree 
models for small m would have been effective. 

Yet another very appealing feature of BART is that it appears to be very 
robust to small changes in the prior and to the choice of m, the number of trees. 
This robustness is illustrated in in Figures 6(a) and (b) which display the in- 
and out-of-sample RMSE obtained by BART as (is, q, k, m) are varied. These 
are based on posterior mean estimates of f(x) from 5000 BART MCMC draws 
(after skipping 1000 burn- in iterations). In each plot of RMSE versus m, the 
plotted text indicates the values of (is, q, k): k = 1, 2 or 3 and (is, q) = d, a or c 
( default /agressive/ conservative). Three striking features of the plot are apparent: 
(i) a very small number of trees (m very small) gives poor results, (ii) as long 
as k > 1, very similar results are obtained from different prior settings, and (iii) 
increasing the number of trees well beyond the number needed to capture the fit, 
results in only a slight degradation of the performance. 

As Figure 6 suggests, the BART fitted values are remarkably stable as the set- 
tings are varied. Indeed, in this example, the correlations between out-of-sample 
fits turn out to be very high, almost always greater than .99. For example, the 
correlation between the fits from the (is, q, k, m)=(3,.9,2,100) setting (a reason- 
able default choice) and the (10, .75, 3, 100) setting (a very conservative choice) 
is .9948. Replicate runs with different seeds are also stable: The correlation 
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Figure 6: BART's robust RMSE performance as (u,q,k,m) is varied: (a) in-sample RMSE 
comparisons and (b) out-of-sample RMSE comparisons. 
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between fits from two runs with the (3, .9, 2, 200) setting is .9994. Such stability 
enables the use of one long MCMC run. In contrast, some models such as neural 
networks require multiple starts to ensure a good optimum has been found. 

5.2.2 Out-of-Sample Comparisons with Competing Methods 

To gauge how well BART performs on the Friedman setup, we compared the out- 
of-sample performance of BART-cv and BART-default with the four competing 
methods considered in Section 5.1, the Lasso, gradient boosting, neural nets and 
random forests. We also considered MARS (Friedman (1991), implemented as 
polymars in R by Kooperberg, Bose & Stone (1997)). 

The methods were compared using 50 independent replications of the simula- 
tion used in Section 5.2.1, namely 100 independent values of (x,y) from (26) and 
(27) with p = 10. Each method was then trained on these 100 in-sample values to 
obtain an estimate / of / using (18). With the exception of BART-default, 5-fold 
cross-validation was used to select from the operational parameter values listed 
in Table 1 or for MARS, to select a value of its gcv parameter from 1,...,8. We 
then simulated an additional 1000 out-of-sample x values from (26). For each of 
our 50 training samples, the performance of each method was then evaluated by 
the RMSE = ^E?=i(/(^i) - /Oi)) 2 over the n = 1000 out-of-sample values. 

Average RMSEs over 50 replicates and standard errors of averages are given 
in Table 3. All the methods explained substantial variation, as the average 
RMSE for the constant model (y = y) is 4.87. However, both BART-cv and 
BART-default substantially outperformed all the other methods by a significant 
amount. Surprisingly, BART-default performed here even better than BART-cv. 
Confirming what we saw in Section 5.1, BART- default's simplicity and speed 
make it an ideal tool for automatic exploratory investigation. Finally, we note 
that BART-cv chose the default (u, q, k) = (3, 0.90, 2.0) most frequently (20% of 
the replicates). 
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Method 



average RMSE se(RMSE) 



Random Forests 

Linear Regression 

Neural Nets 

Boosting 

MARS 

BART-cv 

BART-default 



2.655 0.025 

2.618 0.016 

2.156 0.025 

2.013 0.024 

2.003 0.060 

1.787 0.021 

1.759 0.019 



Table 3: Out-of-sample performance on 50 replicates of the Friedman data. 

5.2.3 Finding Low Dimensional Structure in High Dimensional Data 

Of the p variables Xi,...,x p from (26), / in (27) is a function of only five 
X\,...,X5. Thus the problem we have been considering is one of drawing in- 
ference about a five dimensional signal embedded in a p dimensional space. In 
the previous subsection we saw that when p = 10, the setup used by Friedman 
(1991), BART could easily detect and draw inference about this five dimensional 
signal with just n = 100 observations. We now consider the same problem with 
substantially larger values of p to illustrate the extent to which BART can find 
low dimensional structure in high dimensional data. For this purpose, we re- 
peated the analysis displayed in Figure 3 with p = 20, 100 and 1000 but again 
with only n = 100 observations. We used BART with the same default setting 
of (z/, q, k) = (3, 0.90, 2) and m = 100 with one exception; we used the naive es- 
timate a (the sample standard deviation of Y) rather the least squares estimate 
to anchor the qth prior quantile to allow for data with p > n. Note that because 
the naive a is very likely to be larger than the least squares estimate, it would 
also have been reasonable to use the more aggressive prior setting for (u, q). 

Figure 7 displays the in-sample and out-of-sample BART inferences for the 
larger values p = 20, 100 and 1000. The in-sample estimates and 90% posterior 
intervals for f(x) are remarkably good for every p. As would be expected, the 
out-of-sample plots show that extrapolation outside the data becomes less reliable 
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Figure 7: Inference about Friedman's function in p dimensions. 
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as p increases. Indeed the estimates stray further from the truth especially at the 
boundaries, and the posterior intervals widen (as they should). Where there is less 
information, it makes sense that BART pulls towards the center because the prior 
takes over and the /x's are shrunk towards the center of the y values. Nonetheless, 
when the dimension p is so large compared to the sample size n = 100, it is 
remarkable that the BART inferences are at all reliable, at least in the middle of 
the data. 

In the third column of Figure 7, it is interesting to note what happens to the 
MCMC sequence of a draws. In each of these plots, the solid line at a = 1 is the 
true value and the dashed line at a — 4.87 is the naive estimate used to anchor 
the prior. In each case, the a sequence repeatedly crosses a — 1. However as p 
gets larger, it increasingly tends to stray back towards larger values, a reflection 
of increasing uncertainty. Lastly, note that the sequence of a draws in Figure 7 
are systematically higher than the o draws in Figure 3(c). This may be due in 
part to the fact that the regression a rather than the naive a was used to anchor 
the prior in Figure 3. Indeed if the naive a was instead used for Figure 3, the a 
draws would similarly rise. 

A further attractive feature of BART is that it appears to avoid being misled 
by pure noise. To gauge this, we simulated n = 100 observations from (26) with 
/ = for p = 10, 100, 1000 and ran BART with the same settings as above. With 
p — 10 and p = 100 all intervals for / at both in-sample and out-of-sample x 
values covered or were close to clearly indicating the absence of a relationship. 
At p = 1000 the data becomes so uninformative that our prior, which suggests 
that there is some fit, takes over and some in-sample intervals are far from 0. 
However, the out-of-sample intervals still tend to cover and are very large so 
that BART still indicates no evidence of a relationship between y and x. 



33 



5.3 Classification: A Drug Discovery Application 

Our last example illustrates an application of the BART probit approach de- 
scribed in Section 4 to a drug discovery classification problem. In such problems, 
the goal is to predict the "activity" of a compound against a biological target, 
using predictor variables that characterize the molecular structure of the com- 
pound. By "activity" , one typically means the ability to effect a desired outcome 
against some biological target, such as inhibiting or killing a certain virus. 

The data we consider describe p = 266 molecular characteristics of n = 29, 374 
compounds, 542 of which were classified as active. These predictors represent 
topological aspects of molecular structure. This data set was collected by the 
National Cancer Institute, and is described in Feng, Lurati, Ouyang, Robinson, 
Wang, Yuan & Young (2003). For the in-sample and out-of-sample comparisons 
below, we randomly split the data into nonoverlapping train and test sets, each 
with 14687 compounds of which 271 were active. 

Designating the activity of a compound by a binary variable (Y — 1 if ac- 
tive and Y = otherwise), we applied BART probit to the train set to obtain 
posterior mean estimates of the conditional probabilities, P[Y = l\x], for each 
x vector of the 266 molecular predictor values. For this purpose we used BART 
probit with the default k = 2 and m = 50 trees [y and q have no meaning for 
the probit model, since the latent response is assumed to have variance 1). To 
gauge MCMC convergence, we performed four independent repetitions of 250,000 
MCMC iterations and obtained essentially the same results each time. 

To get a feel for the extent to which BART's posterior mean estimates of 
P[Y = l\x] can be used to identify promising drugs, Figure 8 plots the 20 largest 
P\Y — 1 1 x] estimates for the train and the test sets. Also provided are the 90% 
posterior intervals which convey uncertainty and the identification whether the 
drug was in fact active (y = 1) or not (y = 0). As there are four inactives in each 
plot, the true positive rates in both the train and test sets for these 20 largest 
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Figure 8: BART posterior intervals for the 20 compounds with highest predicted activity, 
using train (a) and test (b) sets. 

estimates is 16/20 = 80%, an impressive gain over the 271/14687 = 1.85% base 
rate. It may be of interest to also note that the test set intervals are slightly 
wider, with an average width of 0.50 compared to 0.47 for the training intervals. 

Going beyond the performance of the top 20 estimates, we assess overall out- 
of-sample predictive performance with Receiver Operating Characteristic (ROC) 
curves based on the orderings of the obtained test set P[Y = 1 1 x] estimates. 
As a frame of reference for comparison, we also applied BART probit with m = 
1 (which we denote by BCART because it is essentially a version of Bayesian 
CART), and a greedy nonBayesian recursive partitioning method implemented as 
the rpart package in R. For RPART, a large tree (around 80 terminal nodes) was 
grown without pruning, as the superior performance of large trees over pruned 
trees on this problem has been previously documented Wang (2005). 

The ROC curves for BART, BCART and RPART, calculated with the R0CR 
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Figure 9: Test set ROC curves for BART, BCART, and RPART. 

package in R and displayed in Figure 9, plot the tradeoffs between true positive 
and false positive rates across thresholds for each method. The AUC (area under 
the curve) values for the BART, BCART and RPART curves here are 0.77, 0.70 
and 0.66, respectively. Noting that higher curves indicate better performance, 
and that a classifier's AUC value is the probability that it will rank a randomly 
chosen y = 1 example higher than a randomly chosen y = example (Bamber 
1975), we see that BART has provided more reliable predictions than BCART 
which in turn has provided more reliable predictions than RPART. We note in 
passing that all three methods are substantially better than random guessing 
which corresponds to a 45-degree line ROC curve with an AUC of 0.50. 

Finally, we turn to the issue of variable selection and demonstrate that by 
decreasing the number of trees m, BART probit can be used, just as BART in 
Section 5.3.1, to identify those predictors which have the most influence on the 
response. For this purpose, we modify the data setup as follows: instead of hold- 
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Figure 10: Variable importance measure, drug discovery example. Values are given for 5, 10 
and 20 trees in the ensemble, for all 266 variables (a) and the 25 variables with the highest 
mean usage (b). Vertical lines in (a) indicate variables whose percent usage exceeds the 95th 
percentile. The 95th percentile is indicated by a horizontal line. 

ing out a test set, all 542 active compounds and a subsample of 542 inactives were 
used to build a model. Four independent chains, each with 1,000,000 iterations, 
were used. The large number of iterations was used to ensure stability in the 
"percent usage" variable selection index (20). Ensembles with 5, 10, and 20 trees 
were considered. 

The results are displayed in Figure 10. The same three variables are selected as 
most important for all three ensembles. Considering that 1/266 ~ 0.004, percent 
usages of 0.050 to 0.100 are quite a bit larger than one would expect if all variables 
were equally important. As expected, variable usage is most concentrated in the 
case of a small ensemble (i.e., m = 5 trees). 
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6 Extensions and Related Work 



Although we have framed BART as a stand alone procedure, it can also be incor- 
porated into larger statistical models for example by adding other components 
such as linear terms or linear random effects. For instance, one might consider a 
model of the form 

Y — hi(x) + h 2 (z) + e, e~iV(0,a 2 ) (28) 

where hi(x) is a sum of trees as in (2) and h 2 {z) is a parametric form involving 
z, a second vector of predictors. One can also extend the sum-of-trees model to 
a multivariate framework such as 

Yi = h i (xi) + € i , (ei,e 2 ,...,e p ) ~iV(0,S), (29) 

where each hi is a sum of trees and E is a p dimensional covariance matrix. If all 
the Xi are the same, we have a generalization of multivariate regression. If the 
Xi are different we have a generalization of Zellner's SUR model (Zellner 1962). 
The modularity of the BART MCMC algorithm in Section 3.1 easily allows for 
such incorporations and extensions. Implementation of linear terms or random 
effects in a BART model would only require a simple additional MCMC step to 
draw the associated parameters. The multivariate version of BART (29) is easily 
fit by drawing each h* given {h*}j^i and E, and then drawing E given all the h*. 

An early version of our work on BART (Chipman, George & McCulloch 2007) 
was published in the proceedings of the conference Advances in Neural Informa- 
tion Processing Systems 2006. Based on this and other preliminary technical 
reports of ours, a variety of extensions and applications of BART have begun to 
appear. Zhang, Shih & Muller (2007) proposed SBART an extension of BART 
obtained by adding a spatial component along the lines of (28). Applied to the 
problem of merging data sets, they found that SBART improved over the conven- 
tional census based method. For the predictive modeling problem of TF-DNA 
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binding in genetics, Liu & Zhou (2007) and Zhou & Liu (2007) considered a 
variety of learning methods, including stepwise linear regression, MARS, neural 
networks, support vector machines, boosting and BART. Concluding that "the 
BART method performed best in all cases", they noted BART's "high predic- 
tive power, its explicit quantification of uncertainty and its interpretability" . By 
keeping track of the per sample inclusion rates, they successfully used BART 
to identify some unusual predictors. Zhang & Hardle (2007) independently dis- 
covered the probit extension of BART, which they call BACT, and applied it 
to credit risk data to predict the insolvency of firms. They found BACT to 
outperform the logit model, CART and support vector machines. Abu-Nimeh, 
Nappa, Wang & Nair (2008) also independently discovered the probit extension 
of BART , which they call CBART, and applied it for the automatic detection 
phishing emails. They found CBART to outperform logistic regression, ran- 
dom forests, support vector machines, CART, neural networks and the original 
BART. Abreveya & McCulloch (2006) applied BART to hockey game penalty 
data and found evidence of referee bias in officiating. Finally, Hill & McCulloch 
(2007) apply BART to a counterfactual causality analysis for assessing the ef- 
fectiveness of a job training program based on observational data. They found 
that BART produced more efficient estimates than propensity score matching, 
propensity-weighted estimators and regression adjustment. Without exception, 
these papers provide further evidence for the remarkable potential of BART. 

7 Discussion 

The essential components of BART are the sum-of-trees model, the regulariza- 
tion prior and the backfitting MCMC algorithm. As opposed to the Bayesian 
approaches of CGM98 and Denison et al. (1998), where a single tree is used to 
explain all the variation in y, each of the trees in BART accounts for only part 
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of the overall fit. This is accomplished with a regularization prior that shrinks 
the tree effects towards a simpler fit. To facilitate the implementation of BART, 
the prior is formulated in terms of rapidly computable forms that are controlled 
by interpretable hyperparameters, and which allow for a highly effective default 
version for immediate "off-the-shelf" use. Posterior calculation is carried out by 
a tailored backfitting MCMC algorithm that appears to converge quickly, effec- 
tively obtaining a (dependent) sample from the posterior distribution over the 
space of sum-of-trees models. A variety of inferential quantities of interest can 
be obtained directly from this sample. 

The application of BART to a wide variety of data sets and a simulation ex- 
periment (Section 5) served to demonstrate many of its appealing features. In 
terms of out-of sample predictive RMSE performance, BART compared favorably 
with boosting, the lasso, MARS neural nets and random forests. In particular, 
the computationally inexpensive and easy to use default version of BART per- 
formed extremely well. In the simulation experiments, BART obtained reliable 
posterior mean and interval estimates of the true regression function as well as 
the marginal predictor effects. BART's performance was seen to be remarkably 
robust to hyperparameter specification, and remained effective when the regres- 
sion function was buried in ever higher dimensional spaces. BART was also seen 
to be a new effective tool for model-free variable selection. Finally, a straight- 
forward probit extension of BART for classification of binary Y was seen to be 
an effective tool for discovering promising drugs on the basis of their molecular 
structure. 
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